# attach packages

library(survival)
library(foreign)
library("survminer")
library('zoo')
library('flexsurvreg')
library('tidyr')
library('eha')
library('dplyr')
library(broom)
library('dotwhisker')
library('simPH')
library(RCurl)

# call datasets. Note that not all are originally in .csv format and will require extraction from zip, then conversion. 

#singh and way -- location: https://journals.sagepub.com/doi/suppl/10.1177/0022002704269655 
sw <- read.csv("singh_and_way/sw.csv")
#Kroenig -- location: https://journals.sagepub.com/doi/suppl/10.1177/0022002708330287 
kro <- read.csv("Kroenig/kroenig.csv")
#Bleek and lorber -- location: https://journals.sagepub.com/doi/suppl/10.1177/0022002713509050 
ble <- read.csv("Bleek/ble.csv")
#Jo and gartzke -- location: http://dss.ucsd.edu/~egartzke/data/jo_gartzke_0207_replicate_0906.dta 
jog <- read.csv("jo_and_gartzke/jog.csv")
#Way and weeks -- location: https://dataverse.harvard.edu/file.xhtml?persistentId=doi:10.7910/DVN/VWBP3L/RYXCBB&version=2.1 
ww <- read.table("way_and_weeks/WayWeeksAJPS.tab", sep = '\t',header = TRUE)
#Reiter -- location: https://dataverse.harvard.edu/file.xhtml?persistentId=doi:10.7910/DVN/HK6PKL/GPIST0&version=1.0 
dr <- read.table("Reiter/reiterfpanucacq.tab", sep = '\t',header = TRUE)
#Montgomery -- location: https://people.reed.edu/~ahm/Projects/ProlifTK/Montgomery2013Stop.zip 
mon <- read.csv("Montgomery/Montgomery 2013.csv")

# analysis

                                                    #@# 1. Singh and Way (2004) #@#

#1a clean - make factors
sw$industry1 <- as.factor(sw$industry1)
sw$allies <- as.factor(sw$allies) 


#1b specify models to check cox ph vs weibull distribution. Model from Singh and Way (2004,873)
#function 'coxreg' cannot be used with clustering: https://r.789695.n4.nabble.com/coxreg-vs-coxph-time-dependent-treatment-td3806161.html 
#first, create risk1 - risk + 1 for survreg
sw$risk1 <- sw$risk+1
sw_weibull_acq <- phreg(Surv(risk1, acquire) ~ cluster(cowcc) +gdpcap + g2 + industry1 + rivalry + centdems + dch5 + rivalry + disputes + openness + chopen5 + allies + polity, dist="weibull", data=sw)
sw_cox_acq1 <- coxreg(Surv(risk1, acquire) ~ gdpcap + g2 + industry1 + rivalry + centdems + dch5 + rivalry + disputes + openness + chopen5 + allies + polity, data=sw)

#1c graphical check on the assumption that Weibull is appropriate. Using the 'coxreg' function is used to generate the model because it permits the check.dist function. Later I use coxph to generate substantive findings. The outputs match.
check.dist(sw_cox_acq1, sw_weibull_acq)
#FINDING 1: parametric model inappropriate, reported in text in TABLE 1.  

#1d Check on proportional hazards
# Documentation: https://r.789695.n4.nabble.com/Why-is-transform-quot-km-quot-the-default-for-cox-zph-td797535.html 
# Important here to have coxph with clustering, greatly impacts the proportional hazards test. 
sw_cox_acq2 <- coxph(Surv(risk1, acquire) ~ gdpcap + g2 + industry1 + centdems + dch5 + rivalry + disputes + openness + chopen5 + allies + polity + cluster(cowcc), robust=T, data=sw)
sw_test <- cox.zph(sw_cox_acq2)

# Sometimes the line above throws an error. This appears to depend on the version of R used with recent versions more likely to cause a problem. Stable version 3.5.2 works OK. Likely this is due to the collinearity extant in the model
# especially related to the economic variables 'g2' 'gdpcap' and 'industry1'. Altering the model to reflect this collinearity (by removing g2 and industry1), gives 'sw_cox_acq3'. I haven't found the cox.zph test for this model to
# cause an error. It also reports violated covariates -- although only 2. This collinearity issue is discussed in the text in final paragraph of section 4. 
#sw_cox_acq3 <- coxph(Surv(risk1, acquire) ~  gdpcap + centdems + dch5 + rivalry + disputes + openness + chopen5 + allies + polity + cluster(cowcc), robust=T, data=sw)
#sw_test2 <- cox.zph(sw_cox_acq3)

#1e  Plot schoenfeld residuals, which assess the proportional hazards by variable
ggcoxzph(sw_test2)
#FINDING 2: according to these tests, 8 variables strongly depend on time. And the model fails the global test (ie the model depends on time). Reported in TABLE 1.

#1f tests for ph in model for pursue. Print and graph results. 
sw_cox_purs <- coxph(Surv(risk1, pursue) ~ gdpcap + g2 + industry1 + centdems + dch5 + rivalry + disputes + openness + chopen5 + allies + polity + cluster(cowcc), robust=T, data=sw)
sw_test_purs <- cox.zph(sw_cox_purs)
sw_test_purs
ggcoxzph(sw_test_purs)
#FINDING 3: tests show that model of pursuit is less affected by time dependence than model of acquisition. Reported in section 4, para 8.

                                                            #@# 2. Kroenig (2009) #@# 

#2a clean data
# make categorical variables factors
# Left censor for existence of program, per Kroenig (2009)
kro$nucass 		<- as.factor(kro$nucass)
kro$industry 	<- as.factor(kro$industry)
kro$rivalry 	<- as.factor(kro$rivalry)
kro$allies 		<- as.factor(kro$allies)
kro$program 	<- as.factor(kro$program)

#2b specify models to check cox ph vs weibull distribution
# Graphical check on the assumption that weibull is appropriate. 
kro_weibull_acq <- phreg(Surv(risk1, acquire) ~ gdpcap + g2 + industry + rivalry + nucass + openness + chopen5 + allies + polity + cluster(cowcc), dist="weibull", data=kro)
kro_cox_acq1 <- coxreg(Surv(risk1, acquire) ~ gdpcap + g2 + industry + rivalry + nucass  + openness + chopen5 + allies + polity + cluster(cowcc), data=kro)
check.dist(kro_cox_acq1, kro_weibull_acq)
#FINDING 4: parametric model does not fit data, reported in TABLE 1.  

#2c. reproduce main model, Kroenig (2009, 174)
kro_cox_acq2 <- coxph(Surv(risk1, acquire) ~ nucass + chopen5 + openness + gdpcap + g2 + industry + rivalry  + allies + polity + cluster(cowcc), data=kro, robust=T)

#2d. Check on proportional hazards assumption. 
kro_phtest <- cox.zph(kro_cox_acq2)
kro_phtest
ggcoxzph(kro_phtest)
#FINDING 5: 5 covariates depend on time, including 'nucass'. The model fails the global test of time dependence. Test graphically by plotting Schoenfeld residuals. Reported in TABLE 1.

#If the collinearity problem reported above for the Singh and Way model persists here, simply run a modified model (drop g2 and industry1) and retest, as above.
#kro_cox_acq2.1 <- coxph(Surv(risk1, acquire) ~ nucass + chopen5 + openness + gdpcap + industry + rivalry  + allies + polity + cluster(cowcc), data=kro, robust=T)
#kro_phtest2 <- cox.zph(kro_cox_acq3)
#kro_phtest2

#2e Generate (log-time) interactions for key variable: nucass. 
cols.integer2 <- c("nucass")
kro[cols.integer2 ] <- sapply(kro[cols.integer2 ],as.integer)
sapply(kro, class)
KrobaseVars <- c('nucass')
kro2 <- tvc(kro, b = KrobaseVars, tvar ='risk1', tfun ='log')

#2f Rerun the model with control for time
kro_cox_acq3 <- coxph(Surv(risk1, acquire) ~ nucass + nucass_log + gdpcap + g2 + industry + rivalry + allies + polity + cluster(cowcc), data=kro2, robust=T, iter.max=5000)

#2g models of nuclear pursuit. 'prog_start' added manually based on 'prog' variable in Kroenig's original dataset. 
kro_cox_purs <- coxph(Surv(risk1, prog_start) ~ cluster(cowcc) + gdpcap + g2 + industry + rivalry + nucass + openness + chopen5 + allies + polity, data=kro, robust=T)
kro_test_purs <- cox.zph(kro_cox_purs)
kro_test_purs
ggcoxzph(kro_test_purs)
#supports FINDING 4: models of pursuit do not have the same time dependence problems as modelling acquisition. Collinearity problem, noted above, re-appears here too.

                                                                #@# 2.  Bleek and Lorber (2010) #@#

#3a clean - Select data from larger ble dataset
# identify factors
# set risk + 1
ble2 <- select(ble, risk1, acquire, pursue, explore, nucsecguarsup, new_econ, nuk7set1, ln_ri_t, nucriv, democ, npt_rati, npt_eff, majpow, regpowt, nucass_full, civass, openness, chopen5, cowcc)
ble2$nucsecguarsup 	<- as.factor(ble2$nucsecguarsup)
ble2$nucriv			<- as.factor(ble2$nucriv)
ble2$npt_rati		<- as.factor(ble2$npt_rati)
ble2$majpow			<- as.factor(ble2$majpow)
ble2$regpowt 		<- as.factor(ble2$regpowt)
ble2$nucass_full	<- as.factor(ble2$nucass_full)
ble2$civass			<- as.factor(ble2$civass)
ble2$risk <- ble2$risk1 + 1 

#3b Reproduce 'core model' acquire. Bleek and Lorber (2013,438). No clustering.
# check coxph vs weibull with graphic
ble_weibull_acq <- phreg(Surv(risk, acquire) ~ nucsecguarsup + new_econ + nuk7set1 + ln_ri_t + nucriv + cluster(cowcc),dist="weibull", data=ble2)
ble_cox_acq1 <- coxreg(Surv(risk, acquire) ~ nucsecguarsup + new_econ + nuk7set1 + ln_ri_t + nucriv + cluster(cowcc), data=ble2)
check.dist(ble_cox_acq1, ble_weibull_acq)
#FINDING 6: parametric model appropriate. TABLE 1.

#3c. Reproduce 'core model' with clustering
ble_cox_acq2 <- coxph(Surv(risk, acquire) ~ nucsecguarsup + new_econ + nuk7set1 + ln_ri_t + nucriv + cluster(cowcc), data=ble2, robust=T)

#3d test proportional hazards assumption, print
ble_phtest1 <- cox.zph(ble_cox_acq2)
ble_phtest1
ggcoxzph(ble_phtest1)
#FINDING 7: tests show 2 covariates including 'nucsecguarsup' depend on time, global test failed. TABLE 1.

#3e Similar check on larger model in Bleek and Lorber (2013,438) confirms suspicions of non-proportional hazard model but result overfitted.
ble_cox_acq3 <- coxph(Surv(risk, acquire) ~ nucsecguarsup + new_econ + nuk7set1 + ln_ri_t + nucriv + nucass_full + civass + npt_rati + npt_eff + majpow + regpowt + democ + openness + chopen5 + cluster(cowcc), data=ble2, robust=T, iter.max=5000)
ble_phtest2 <- cox.zph(ble_cox_acq3)
ble_phtest2

# 3f Generate (log-time) interactions for relevant covariate - nucsecguarsup 
cols.integer2 <- c("nucsecguarsup")
ble2[cols.integer2 ] <- sapply(ble2[cols.integer2 ],as.integer)
sapply(ble2, class)
BlebaseVars <- c('nucsecguarsup')
ble3 <- tvc(ble2, b = BlebaseVars, tvar ='risk1', tfun ='log')

#3g rerun the model with controls
ble_cox_acq4 <- coxph(Surv(risk, acquire) ~ nucsecguarsup + ln_t_nucsecguarsup + new_econ + nucriv + ln_ri_t + nuk7set1 + cluster(cowcc), data=ble3, robust=T, iter.max=5000)
ble_cox_acq4

#3j. models of nuclear assistance, using Bleek's dataset - not reported on in text.
ble_cox_acq5 <- coxph(Surv(risk1, acquire) ~ gdpcap + g2 + industry1 + rivalry + NCAtodate  + allies + polity + cluster(cowcc), data=ble, robust=T)
ble_phtest2 <- cox.zph(ble_cox_acq5)
ble_phtest2
ble$ln_t_nca 	<- as.numeric(ble$NCAtodate)	* log(ble$risk)
ble_cox_acq5 <- coxph(Surv(risk1, acquire) ~ gdpcap + g2 + industry1 + rivalry + NCAtodate + ln_t_nca + allies + polity + cluster(cowcc), data=ble, robust=T)

#bleek model of pursuit. 
ble_cox_purs <- coxph(Surv(risk1, pursue) ~ nucsecguarsup + new_econ + nuk7set1 + ln_ri_t + nucriv + cluster(cowcc), data=ble2, robust=T)
ble_test_purs <- cox.zph(ble_cox_purs)
ggcoxzph(ble_test_purs)
ble_test_purs
#supports FINDING 4 in section 4, para 8 that models of pursuit less affected by time dependence

                                                            #@# 4. Jo and Gartzke (2009) #@#

#4a clean - add risk, factors, censor
jog$risk1 <- jog$year - 1939
jog2 <- jog[jog$nuk_apl==1, ]
jog2$r_nukep	<- as.factor(jog2$r_nukep)
jog2$nuke_a_d	<- as.factor(jog2$nuke_a_d)
jog2$majpow	    <- as.factor(jog2$majpow)
jog2$regpowt	<- as.factor(jog2$regpowt)
jog2$npt_rati 	<- as.factor(jog2$npt_rati)
#Note that 'exit' is used after being inserted manually

#4b Model jog_weibull_acq1. But fails to converge -- problem seems to be with the 'd_isol' and 'ln_xst1' which are highly correlated with other covariates. These are non-significant in Jo and Gartzke, have little basis in theory and are removed by Bleek (2010), I do the same.  
jog_weibull_acq1 <- phreg(Surv(risk1, nuke_exit) ~ nuk7set1 + new_econ + ln_ri_t + nuke_a_d + d_isol + ln_xst1 + democ + npt_eff + majpow + cluster(ccode1),dist="weibull", data=jog2)
jog_model_frame <- select(jog2,nuk7set1 , new_econ , ln_ri_t , nuke_a_d , d_isol , ln_xst1 , democ , npt_eff , majpow)
nums <- unlist(lapply(jog_model_frame, is.numeric))  
jog2_set <- jog_model_frame[, nums]
matriz_cor <- cor(jog2_set)

#4c Model minus 'd_isol' and 'ln_xst1'. Weibull vs coxph. Graphical check.
jog_weibull_acq2 <- phreg(Surv(risk1, nuke_exit) ~ nuk7set1 + new_econ + ln_ri_t + nuke_a_d + democ + npt_eff + majpow + cluster(ccode1),dist="weibull", data=jog2)
jog_cox_acq1 <- coxreg(Surv(risk1, nuke_exit) ~ nuk7set1 + new_econ + ln_ri_t + nuke_a_d + democ + npt_eff  + majpow + cluster(ccode1),  data=jog2)
check.dist(jog_cox_acq1, jog_weibull_acq2)
# FINDING 8: parametric inappropriate. TABLE 1. 

#4d Model including cowcc clustering
jog_cox_acq2 <- coxph(Surv(risk1, nuke_exit)~ nuk7set1 + new_econ + ln_ri_t + nuke_a_d + democ + npt_eff  + majpow + cluster(ccode1), data=jog2, robust=T)

#4e Proportional hazards check
jog_phtest1 <- cox.zph(jog_cox_acq2)
jog_phtest1
ggcoxzph(jog_phtest1)
# FINDING 9: according to the tests, 1 variable depends on time, global test failed. Show graphically.  


#4f test model of pursuit 
jog_cox_purs <- coxph(Surv(risk1, jog$nuk_a_p2)~ nuk7set1 + new_econ  + ln_ri_t +  nuke_a_d + democ + npt_eff  + majpow + cluster(ccode1), data=jog, robust=T)
jog_test_purs <- cox.zph(jog_cox_purs)
ggcoxzph(jog_test_purs)
jog_test_purs
#Less support for FINDING 4 since 1 variable violates ph, and global test failed. 

                                                                    #@#5. Way and Weeks (2014) (pursue only, not covered in text) #@#

#5a clean - factors and risk
ww$persdumjlw_lag <- as.factor(ww$persdumjlw_lag)
ww$risk1 <- ww$risk+1

# 5b. weibull vs coxph and graphical check
ww_weibull_acq <- phreg(Surv(timesw2, failsw) ~ cluster(ccode) +persdumjlw_lag + lpop + land, dist="weibull", data=ww)
ww_cox_acq1 <- coxreg(Surv(timesw2, failsw) ~ cluster(ccode) +persdumjlw_lag + lpop + land, data=ww)

check.dist(ww_cox_acq1, ww_weibull_acq)

# 5c. Build model with clustering and strata, which is used by Way and Weeks in their paper. 
ww_cox_acq2 <- coxph(Surv(timesw2, failsw) ~ persdumjlw_lag + lpop + land + cluster(stratsw), data=ww)

# 5d.Test proportional hazard
ww_test <- cox.zph(ww_cox_acq2)

# 5e. graph Schoenfeld residuals 
ggcoxzph(ww_test)
#Support for FINDING 4

                                                                  #@# 6. Reiter (2014) #@#

#6a. clean -  factors
dr$industry		<- as.factor(dr$industry)
dr$rivalry	 	<- as.factor(dr$rivalry)
dr$usnucdp		<- as.factor(dr$usnucdp)
dr$nucass		<- as.factor(dr$nucass)

#6b. cox vs weibull. graphical check
dr_weibull_acq<- phreg(Surv(risk1, acquire) ~ nucass + gdpcap + g2 + industry + polity + usnucdp + usnuc + ustroops + cluster(cowcc),dist="weibull", data=dr)
dr_cox_acq1 <- coxreg(Surv(risk1, acquire) ~ nucass + gdpcap + g2 + industry + polity + usnucdp + usnuc + ustroops + cluster(cowcc), data=dr)
check.dist(dr_cox_acq1, dr_weibull_acq)
#6c. FINDING 10: parametric model appropriate. TABLE 1.

#6c. Model with clustering
dr_cox_acq2 <- coxph(Surv(risk1, acquire)~ nucass + gdpcap + g2 + industry + polity + usnucdp + usnuc + ustroops + cluster(cowcc), data=dr, robust=T)

#6d. check Proportional hazards. graphical test. 
dr_phtest1 <- cox.zph(dr_cox_acq2)
dr_phtest1
ggcoxzph(dr_phtest1)
# FINDING 11: 4 covariates violate ph, including 'usnucdp', global test failed. TABLE 1. 

#6e. generate (log-time) interactions for relevant covariates (usnucdp1, usnuc)
dr$ln_t_usnucdp1		<- as.numeric(dr$usnucdp)			* log(dr$risk1)
dr$ln_t_usnuc			<- as.numeric(dr$usnuc)				* log(dr$risk1)

# 6f. Model 'dr_cox_acq3' sometimes fails to converge - likely due to the 'usnuc' variable which is closely correlated with aquire, drop variable and rerun. Note that effect of key variable usnucdp does not change substantively between model 3 and 4.
# Another collinearity issue similar to that mentioned in section 4, para 8
dr_cox_acq3 <- coxph(Surv(risk1, acquire)~ nucass+ gdpcap + g2 + industry + polity  + usnucdp + ln_t_usnucdp1 + usnuc + ln_t_usnuc +  ustroops + cluster(cowcc), data=dr, robust=T)
dr_cox_acq3
dr_cox_acq4 <- coxph(Surv(risk1, acquire)~ nucass + gdpcap + g2 + industry + polity  + usnucdp + ln_t_usnucdp1 + cluster(cowcc), data=dr, robust=T)
dr_cox_acq4

                                                                    #@# 7. Montgomery (2013) #@# 

# 7a. clean - factors and risk.
mon$industry1	<- as.factor(mon$industry1)
mon$nucass		<- as.factor(mon$nucass)
mon2 <- mon[(mon$risk3==0), ]


#7b. weibull vs coxph. graphical check. 
mon_weibull_acq<- phreg(Surv(risk3, level3onset) ~ nucass + neopat1 + neoass + industry1 + cluster(cowcc),dist="weibull", data=mon)
mon_cox_acq1 <- coxreg(Surv(risk3, level3onset) ~ nucass + neopat1 + neoass + industry1, data=mon)
check.dist(mon_cox_acq1, mon_weibull_acq)
# FINDING 12: parametric inappropriate. TABLE 1.

# 7c. Reproduce model A7 with clustering 
mon_cox_acq2 <- coxph(Surv(risk3, level3onset) ~ nucass + neopat1 + neoass + industry1 + cluster(cowcc), data=mon, robust = T)
mon_cox_acq2

#7d. test for proportional hazards. Check graphically. 
mon_phtest1 <- cox.zph(mon_cox_acq2)
mon_phtest1
ggcoxzph(mon_phtest1)
#FINDING 13: 1 variable violates ph. TABLE 1. 

#7e Generate (log-time) interactions for relevant covariate - neopatrimonial
mon$ln_t_neopat1		<- as.numeric(mon$neopat1)			* log(mon$risk1+1)
mon_cox_acq3 <- coxph(Surv(risk1, level3onset) ~neopat1*NCAtodate + ln_t_neopat1 + cluster(cowcc), data=mon, robust = T)
mon_cox_acq3


#8 PLOTS

#8a Simulate hazard ratios for Kroenig, Montgomery's version of Fuhrmann NCA test, Bleek, Reiter

Simkro2 <- coxsimtvc(kro_cox_acq3, b = "nucass", btvc = "nucass_log", tfun = "log", from = 5, to = 56, Xj = 1, by = 1, ci = 0.9, nsim = 100, qi = "Hazard Ratio")
Simmon2 <- coxsimLinear(mon_cox_acq3, b = "NCAtodate", Xj = 1:0, nsim = 100,ci=0.9, qi= 'Hazard Ratio')
Simble2 <- coxsimtvc(ble_cox_acq4, b = "nucsecguarsup", btvc = "ln_t_nucsecguarsup", tfun = "log", from = 1, to = 56, Xj = 1, by = 1, ci = 0.9, nsim = 100, qi = "Hazard Ratio")
Simdr2 <- coxsimtvc(dr_cox_acq4, b = "usnucdp1", btvc = "ln_t_usnucdp1", tfun = "log", from = 10, to = 56, Xj = 1, by = 1, ci = 0.9, nsim = 100, qi = "Hazard Ratio")

#8b Generate and save plots for kroenig and bleek findings 

figure1 <- simGG(Simkro2, spalette = 'Set2', xlab = 'Time (years)', ylab = 'Hazard ratio of acquiring nuclear weapon', legend = F) + ylim(0,2000)
figure1
ggsave("figure1.jpeg", dpi = 300)

figure2 <- simGG(Simble2, spalette = 'Set2', xlab = 'Time (years)', ylab = 'Hazard ratio of acquiring nuclear weapon', legend = F) + ylim(0,10)
figure2
ggsave("figure2.jpeg", dpi = 300)

#8c. FINDINGs 14&15: nucass and nucsecguarall depend on time. FIGURES 1 and 2.


#8c Simulate hazard ratios for original models

Simkro1 <- coxsimLinear(kro_cox_acq2, b = "nucass", Xj = 1:0, nsim = 100, ci=0.9, qi= 'Hazard Ratio')
Simmon1 <- coxsimLinear(mon_cox_acq3, b = "NCAtodate", Xj = 0:20, nsim = 100, ci=0.9, qi= 'Hazard Ratio')
Simble1 <- coxsimLinear(ble_cox_acq2, b = "nucsecguarsup1", Xj = 1:0, nsim = 100, ci=0.9, qi= 'Hazard Ratio')
Simdr1 <- coxsimLinear(dr_cox_acq2, b = "usnucdp1", Xj = 1:0, nsim = 100, ci=0.9, qi= 'Hazard Ratio')

Simkro1 <- Simkro1[!(Simkro1$QI >1000),]
Simkro2 <- Simkro2[!(Simkro2$QI >1000),]
Simmon1 <- Simmon1[!(Simmon1$QI >1000),]
Simmon2 <- Simmon2[!(Simmon2$QI >1000),]
Simble1 <- Simble1[!(Simble1$QI >1000),]
Simble2 <- Simble2[!(Simble2$QI >1000),]
Simdr1 <- Simdr1[!(Simdr1$QI >1000),]
Simdr2 <- Simdr2[!(Simdr2$QI >1000),]

#8d Get mean simulated hazard ratios 

kromean1 <- mean(Simkro1$QI)
kromedian2 <- mean(Simkro2$QI)
monmean1 <-  mean(Simmon1$QI)
monmean2 <- mean(Simmon2$QI)
blemean1 <- mean(Simble1$QI)
blemedian2<- mean(Simble2$QI)
drmean1<- mean(Simdr1$QI)
drmedian2<- mean(Simdr2$QI)

#8e Generate plots for kroenig, montgomery, bleek, reiter

simulationskro <- tibble(model = c('Kroenig', 'Kroenig (new)'), author = c('Kroenig', 'Kroenig'), term = c('   Kroenig (2009)', 'time adjusted'),
                         estimate=c(kromean1, kromedian2),
                         conf.high=c(quantile(Simkro1$QI,0.05), quantile(Simkro2$QI,0.95)),
                         conf.low=c(quantile(Simkro1$QI,0.05), quantile(Simkro2$QI,0.05)))
k <- dwplot(simulationskro, dots_arg=list(size=100)) + theme(panel.background = element_rect(fill = "white", colour = "darkgrey"),text = element_text(size=15),
                                                             axis.text.x = element_text(angle=0, hjust=1)) + scale_colour_grey() + theme(legend.position = "none") +xlab("") + geom_vline(xintercept = 1, colour = "darkgrey", linetype = 2) 

simulationsmon <- tibble(model = c('Fuhrmann', 'Fuhrmann (new)'), author = c('Fuhrmann', 'Fuhrmann'), term = c('Montgomery (2009)', 'time adjusted'),
                         estimate=c(monmean1, monmean2),
                         conf.high=c(quantile(Simmon1$QI,0.95), quantile(Simmon2$QI,0.95)),
                         conf.low=c(quantile(Simmon1$QI,0.05), quantile(Simmon2$QI,0.05)))
m <- dwplot(simulationsmon, dots_arg=list(size=100)) + theme(panel.background = element_rect(fill = "white", colour = "darkgrey"),text = element_text(size=15),
                                                             axis.text.x = element_text(angle=0, hjust=1)) + scale_colour_grey() + theme(legend.position = "none") +xlab("") + geom_vline(xintercept = 1, colour = "darkgrey", linetype = 2) 

simulationsble <- tibble(model = c( 'Bleek', 'Bleek (new)'), author = c('Bleek', 'Bleek'), term = c('     Bleek & Lorber (2010)', 'time adjusted'),
                         estimate=c(blemean1, blemedian2),
                         conf.high=c(quantile(Simble1$QI,0.95), quantile(Simble2$QI,0.95)),
                         conf.low=c(quantile(Simble1$QI,0.05), quantile(Simble2$QI,0.05)))
b <- dwplot(simulationsble, dots_arg=list(size=100)) + theme(panel.background = element_rect(fill = "white", colour = "darkgrey"),text = element_text(size=15),
                                                             axis.text.x = element_text(angle=0, hjust=1)) + scale_colour_grey() + theme(legend.position = "none") +xlab("") + geom_vline(xintercept = 1, colour = "darkgrey", linetype = 2) 

simulationsdr <- tibble(model = c( 'reiter', 'reiter (new)'), author = c('reiter', 'reiter'), term = c('   Reiter (2010)', 'time adjusted'),
                        estimate=c(drmean1, drmedian2),
                        conf.high=c(quantile(Simdr1$QI,0.95), quantile(Simdr2$QI,0.95)),
                        conf.low=c(quantile(Simdr1$QI,0.05), quantile(Simdr2$QI,0.05)))
d <- dwplot(simulationsdr, dots_arg=list(size=100)) + theme(panel.background = element_rect(fill = "white", colour = "darkgrey"),text = element_text(size=15),
                                                            axis.text.x = element_text(angle=0, hjust=1)) + scale_colour_grey() + theme(legend.position = "none") +xlab("Hazard ratios") + geom_vline(xintercept = 1, colour = "darkgrey", linetype = 2) 
#8f make and print joint plot 

figure3 <- ggarrange(k, m, b, d,
                    labels = c("1 Sensitive", "2 Civilian", '3 Guarantee', '4 US ally'),
                    ncol = 1, nrow = 4)

figure3
ggsave("figure3.jpeg", dpi = 300)

